#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Feb 17 11:44:57 2024

@author: jcfq2
"""

import os

jwst_dir='/Users/jcfq2/data/observations/jwst'

os.chdir(jwst_dir)
import matplotlib.colors as colours
import matplotlib.pyplot as plt
# from astropy.visualization import make_lupton_rgb
# from astropy.table import Table
# import ch4_fiddlesticks as ch4
# import pandas as pd
# import h3ppy
import glob
# import spiceypy as spice
# from astropy.io import fits
import numpy as np
from importlib import reload
import JWSTSolarSystemPointing as jssp
reload(jssp)
#import jwst_uranus as jwstu
from scipy.interpolate import interp1d

# In[0]

kernel_dir = '/Users/jcfq2/data/observations/jwst/kernels/'
jssp.load_kernels(kdir=kernel_dir)

# # Set up a h3ppy object too, always useful
# h3p_model = h3ppy.h3p()

# In[0] set up new colormaps



colors2 = plt.cm.Greens_r(np.linspace(0, 1, 220))
colors1 = plt.cm.cubehelix_r(np.linspace(0.0, 0.4, 32))

# combine them and build a new colormap
colors = np.vstack((colors1, colors2))

for i in range(32):
    colors[31-i,0]=colors[32,0]*(32-(i+1))/32
    colors[31-i,1]=colors[32,1]*(32-(i+1))/32
    colors[31-i,2]=colors[32,2]*(32-(i+1))/32
# for i in range(32): print(colors[32,1]*(32-(i+1))/32)

mygreen = colours.LinearSegmentedColormap.from_list('my_green', colors)


colors2 = plt.cm.Reds_r(np.linspace(0, 1, 220))
colors1 = plt.cm.cubehelix_r(np.linspace(0.0, 0.4, 32))

# combine them and build a new colormap
colors = np.vstack((colors1, colors2))

for i in range(32):
    colors[31-i,0]=colors[32,0]*(32-(i+1))/32
    colors[31-i,1]=colors[32,1]*(32-(i+1))/32
    colors[31-i,2]=colors[32,2]*(32-(i+1))/32
# for i in range(32): print(colors[32,1]*(32-(i+1))/32)

myred = colours.LinearSegmentedColormap.from_list('my_red', colors)


colors2 = plt.cm.Blues_r(np.linspace(0, 1, 220))
colors1 = plt.cm.cubehelix_r(np.linspace(0.0, 0.4, 32))

# combine them and build a new colormap
colors = np.vstack((colors1, colors2))

for i in range(32):
    colors[31-i,0]=colors[32,0]*(32-(i+1))/32
    colors[31-i,1]=colors[32,1]*(32-(i+1))/32
    colors[31-i,2]=colors[32,2]*(32-(i+1))/32
# for i in range(32): print(colors[32,1]*(32-(i+1))/32)

myblue = colours.LinearSegmentedColormap.from_list('my_blue', colors)



colors2 = plt.cm.Purples_r(np.linspace(0, 1, 220))
colors1 = plt.cm.cubehelix_r(np.linspace(0.0, 0.4, 32))

# combine them and build a new colormap
colors = np.vstack((colors1, colors2))

for i in range(32):
    colors[31-i,0]=colors[32,0]*(32-(i+1))/32
    colors[31-i,1]=colors[32,1]*(32-(i+1))/32
    colors[31-i,2]=colors[32,2]*(32-(i+1))/32
# for i in range(32): print(colors[32,1]*(32-(i+1))/32)

mypurple = colours.LinearSegmentedColormap.from_list('my_purple', colors)







# In[1]


import cartopy.crs as ccrs
crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))

 
files = sorted(glob.glob(
    "/Users/jcfq2/data/observations/jwst/5308_dither_separated/*.fits"))
geo = jssp.JWSTSolarSystemPointing(files[0])
wave = geo.get_wavelength()
#  just to get geo for conversion

# nam=np.load('saturn_spectra_map.npy')
# nac=np.load('saturn_spectra_count.npy')
nam=np.load('saturn_spectra_map_2_los.npy')
nac=np.load('saturn_spectra_count_2_los.npy')
# nam=np.load('saturn_spectra_map_2_narrow.npy')
# nac=np.load('saturn_spectra_count_2_narrow.npy')

scale=2

nam[-90*scale:-45*scale,:,:]=nam[-90*scale:-45*scale,:,:]+nam[0:45*scale,:,:]
nac[-90*scale:-45*scale,:,:]=nac[-90*scale:-45*scale,:,:]+nac[0:45*scale,:,:]
nam[45*scale:90*scale,:,:]=nam[45*scale:90*scale,:,:]+nam[-45*scale:,:,:]
nac[45*scale:90*scale,:,:]=nac[45*scale:90*scale,:,:]+nac[-45*scale:,:,:]



c = 2.99792458e+8

for ww in range(len(wave)):
    nam[:,:,ww]=nam[:,:,ww] / (wave[ww] * 1e-6)**2 * c * 1.0e-26
# spec = spec.copy() / (wave * 1e-6)**2 * c * 1.0e-26

nam=nam*1000 # mW

wavemin_F323 = 3.148
wavemax_F323 = 3.38
whw_F323 = np.argwhere((wave > wavemin_F323) & (wave < wavemax_F323)).flatten()


wavemin_F323f = 3.148
wavemax_F323f = 3.35
whw_F323f = np.argwhere((wave > wavemin_F323f) & (wave < wavemax_F323f)).flatten()


# could be problematic - lots going on in the background 
file_F323='/Users/jcfq2/data/observations/jwst/saturn/PSG_files/psg_rad_F323.txt'
bck_F323 = np.loadtxt(file_F323)
f_F323 = interp1d(bck_F323[:,0], bck_F323[:,1], kind='linear')  # You can change 'linear' to other interpolation methods like 'cubic'
bck_F323_whw = f_F323(wave[whw_F323])

F323_values=np.loadtxt('F323N_filter.txt',delimiter=' ')

f_F323f = interp1d(F323_values[:,0], F323_values[:,1], kind='linear')  # You can change 'linear' to other interpolation methods like 'cubic'
filter_F323_sub = f_F323f(wave[whw_F323f])


filter_F323 = np.zeros_like(wave)
filter_F323[whw_F323f]=filter_F323_sub    
# filter_F323_h3p = filter_F323*1.0
# filter_F323_without_h3p = filter_F323*1.0


bck_F323_whw=filter_F323[whw_F323]*bck_F323_whw







sat_T=nam[:,:,600]/nac[:,:,600]


ab=nam[:,:,2239]+nam[:,:,1568]+nam[:,:,1642]+nam[:,:,1010]+nam[:,:,1011]+nam[:,:,1018]+nam[:,:,1019]
cd=(nam[:,:,2241]*0.8*0.5+nam[:,:,2236]*1.19*0.5)+((nam[:,:,1571]+nam[:,:,1565])*0.5)+(0.4*nam[:,:,1640]/5.21)+(nam[:,:,1004]+nam[:,:,1005]*2)


img_h3p=(ab-cd)/nac[:,:,600]


F323_full_image = np.zeros_like(img_h3p)

for xx in range(img_h3p[:,0].size):
    for yy in range(img_h3p[0,:].size):
        if np.isfinite(nam[xx, yy,600]):
            # print(xx,yy)
            pixel_F323=np.nansum(nam[xx, yy,:]/nac[xx, yy,:]*filter_F323)
            # pixel_F323_h3p=np.nansum(geoim[:,xx,yy]*filter_F323_h3p)
            # pixel_F323_noh3p=np.nansum(geoim[:,xx,yy]*filter_F323_without_h3p)
        
            # F323_h3p_image[xx,yy] = pixel_F323_h3p
            # noF323_h3p_image[xx,yy] = pixel_F323_noh3p
            F323_full_image[xx,yy] = pixel_F323

# F323_full_image=np.rot90(F323_full_image,3)    
    

# Set a wavelength range where the non-LTE CH4 is
whw = np.argwhere((wave > 3) & (wave < 4.6)).flatten()
print(whw)
# Fit the non-LTE CH4 background
# crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))
crs = ccrs.NorthPolarStereo(globe=ccrs.Globe(flattening=(0.0)))

    
aaa=np.nanmean(nam[:,:,672:677]/nac[:,:,672:677],axis=2)
bbb=np.nanmean(nam[:,:,672+7:677+7]/nac[:,:,672+7:677+7],axis=2)
ccc=np.nanmean(nam[:,:,672+22:677+22]/nac[:,:,672+22:677+22],axis=2)
# ddd=np.nanmedian(geoim[whw_331,:,:],axis=0)

img_methane_fun = bbb-aaa
img_methane_hot = ccc-aaa

wavemin = 5.0
wavemax = 5.2
whw = np.argwhere((wave > wavemin) & (wave < wavemax)).flatten()
img_heat = np.nanmean(nam[:, :,whw]/nac[:, :,whw], axis=2)

wavemin = 3.0
wavemax = 3.2
whw = np.argwhere((wave > wavemin) & (wave < wavemax)).flatten()
img_reflect = np.nanmean(nam[:, :,whw]/nac[:, :,whw], axis=2)


# wavemin = 3.21
# wavemax = 3.25
# whw = np.argwhere((wave > wavemin) & (wave < wavemax)).flatten()
# img_323_nofilter = np.nanmean(nam[:, :,whw]/nac[:, :,whw], axis=2)



# h3p = h3ppy.h3p()
# h3p.set(T=500, N=6e14, wave=wave, R=2000)
# h3p_model=h3p.model()
# h3p_model=h3p_model/np.max(h3p_model)

# filter_F323_h3p[h3p_model < 0.001]=0
# filter_F323_without_h3p[h3p_model > 0.001]=0

    

# F323_h3p_image = np.zeros_like(img_reflect)
# noF323_h3p_image = np.zeros_like(ef)


# %%
maxint=np.nanmax(img_h3p[45*scale:405*scale,:].T)
maxint=1

map_h3p=np.fliplr(img_h3p[45*scale:405*scale,:].T)/maxint
diff_h3p=map_h3p+0.
median_h3p=np.nanmedian(map_h3p,axis=1)
for i in range(360*scale): diff_h3p[:,i]=map_h3p[:,i]-median_h3p

map_reflect= np.fliplr(img_reflect[45*scale:405*scale,:].T)/maxint
diff_reflect=map_reflect+0.
median_reflect=np.nanmedian(map_reflect,axis=1)
for i in range(360*scale): diff_reflect[:,i]=map_reflect[:,i]-median_reflect

map_heat= np.fliplr(img_heat[45*scale:405*scale,:].T)/maxint
diff_heat=map_heat+0.
median_heat=np.nanmedian(map_heat,axis=1)
for i in range(360*scale): diff_heat[:,i]=map_heat[:,i]-median_heat

map_methane_fun= np.fliplr(img_methane_fun[45*scale:405*scale,:].T)/maxint
diff_methane_fun=map_methane_fun+0.
median_methane_fun=np.nanmedian(map_methane_fun,axis=1)
for i in range(360*scale): diff_methane_fun[:,i]=map_methane_fun[:,i]-median_methane_fun

map_methane_hot= np.fliplr(img_methane_hot[45*scale:405*scale,:].T)/maxint
diff_methane_hot=map_methane_hot+0.
median_methane_hot=np.nanmedian(map_methane_hot,axis=1)
for i in range(360*scale): diff_methane_hot[:,i]=map_methane_hot[:,i]-median_methane_hot

map_323= np.fliplr(F323_full_image[45*scale:405*scale,:].T)/maxint
diff_323=map_323+0.
median_323=np.nanmedian(map_323,axis=1)
for i in range(360*scale): diff_323[:,i]=map_323[:,i]-median_323




# %%

import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

fig2 = plt.figure(figsize=(12,8),dpi=300)


fig2.subplots_adjust(hspace=0,wspace=0.08)
bbox_props1= dict(boxstyle="circle,pad=0.15", fc="whitesmoke", ec="silver", lw=2)


ax = plt.subplot(2,4,1,projection=ccrs.NorthPolarStereo(central_longitude=180))
ax.set_extent([0, 360, 68, 90], ccrs.PlateCarree())
ish3p=ax.imshow(map_h3p, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap=mygreen,vmin=0.0)
# ax.gridlines()
ax.text(360-315,65,'a', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)
cbar=fig2.colorbar(ish3p,location='bottom',aspect=8,pad=0.01,shrink=0.85)
cbar.ax.set_xlabel('H$_3^+$ [mW/m$^2$/sr]')

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,70,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)

# ax.text(310,40,'a')
# for aa in range(0,310,10): ax.text(aa/100.,aa/100.,'a')

ax = plt.subplot(2,4,2,projection=ccrs.NorthPolarStereo(central_longitude=180))
ax.set_extent([0, 360, 68, 90], ccrs.PlateCarree())
isdh3p=ax.imshow(diff_h3p, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='seismic',vmax=np.nanmax(diff_h3p),vmin=(-1)*np.nanmax(diff_h3p))
# ax.gridlines()
ax.text(360-315,65,'b', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)
for drift in range(4): ax.plot([35+90*drift,35+90*drift],[90,40], transform=ccrs.PlateCarree(),c='hotpink',linestyle='dotted')
cbar=fig2.colorbar(isdh3p,location='bottom',aspect=8,pad=0.01,shrink=0.85)
cbar.ax.set_xlabel('$\Delta$ H$_3^+$ [mW/m$^2$/sr]')

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,70,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)

ax.text(360-230,67,'0$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='mediumvioletred')
ax.text(360-50,67,'180$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="right", va="center",size=8,c='mediumvioletred')
ax.text(360-329,67.5,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='mediumvioletred')
# ax.text(360-330,69,'90$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='hotpink')
ax.text(360-144,67.7,'270$^{\circ}\Psi_N$', transform=ccrs.PlateCarree(), ha="left", va="center",size=8,c='mediumvioletred')



ax = plt.subplot(2,4,3,projection=ccrs.NorthPolarStereo(central_longitude=180))
ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
ish3p2=ax.imshow(map_h3p, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),norm=colours.PowerNorm(1,vmin=1/20,vmax=0.2),cmap="Greens_r")
# ax.gridlines()
ax.text(360-315,30,'c', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)
cbar=fig2.colorbar(ish3p2,location='bottom',aspect=8,pad=0.01,shrink=0.85)
cbar.ax.set_xlabel('H$_3^+$ [mW/m$^2$/sr]')

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,40,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (40,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)


ax = plt.subplot(2,4,4,projection=ccrs.NorthPolarStereo(central_longitude=180))

# ax = plt.subplot(2,4,4,projection=ccrs.NorthPolarStereo(central_longitude=180))
# ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
# isdh3p=ax.imshow(diff_h3p, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='seismic',vmax=0.2*np.nanmax(diff_h3p),vmin=(-0.2)*np.nanmax(diff_h3p))
# ax.gridlines()
# ax.text(360-315,30,'b-ii', transform=ccrs.PlateCarree())
# for drift in range(4): ax.plot([35+90*drift,35+90*drift],[90,40], transform=ccrs.PlateCarree())
# cbar=fig2.colorbar(isdh3p,location='bottom',aspect=8,pad=0.01,shrink=0.85)
# cbar.ax.set_xlabel('$diff H_3^+$')


ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
is323=ax.imshow(map_heat, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),vmin=0,cmap=myred)
# ax.gridlines()

# for drift in range(6): ax.plot([50+60*drift,50+60*drift],[90,40], transform=ccrs.PlateCarree())
ax.text(360-315,30,'d', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)


cbar=fig2.colorbar(is323,location='bottom',aspect=8,pad=0.01,shrink=0.85)
cbar.ax.set_xlabel('5.0-5.2$\mu$m [mW/m$^2$/sr]')

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,40,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (40,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)




ax = plt.subplot(2,4,5,projection=ccrs.NorthPolarStereo(central_longitude=180))


ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
isreflect=ax.imshow(map_reflect, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='Greys_r')
# ax.gridlines()
cbar=fig2.colorbar(isreflect,location='bottom',aspect=8,pad=0.01,shrink=0.85)
cbar.ax.set_xlabel('3.0-3.2$\mu$m [mW/m$^2$/sr]')
ax.text(360-315,30,'e', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,40,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (40,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)


def func(arr1d):
    kernel = np.ones(2*scale) / 2*scale
    return np.convolve(arr1d, kernel, mode='same')

first_pass = np.apply_along_axis(func, axis=0, arr=diff_reflect)
diff_reflect_smooth = np.apply_along_axis(func, axis=1, arr=first_pass)

ax = plt.subplot(2,4,6,projection=ccrs.NorthPolarStereo(central_longitude=180))
ax.set_extent([0, 360, 68, 90], ccrs.PlateCarree())
isdreflect=ax.imshow(diff_reflect_smooth/16, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='seismic',vmin=np.nanmax(diff_reflect)/(-8),vmax=np.nanmax(diff_reflect)/8)
# ax.gridlines()

for drift in range(6): ax.plot([50+60*drift,50+60*drift],[90,40], transform=ccrs.PlateCarree(),c='lawngreen',linestyle='--')
cbar=fig2.colorbar(isdreflect,location='bottom',aspect=8,pad=0.01,shrink=0.85)
cbar.ax.set_xlabel('$\Delta$ 3.0-3.2$\mu$m [mW/m$^2$/sr]')
ax.text(360-315,65,'f', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text(((360-longit+180) % 360),70,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)



ax.plot([100,100],[0,30],color='b',linewidth=5)

ax = plt.subplot(2,4,7,projection=ccrs.NorthPolarStereo(central_longitude=180))


ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
isfun=ax.imshow(map_methane_fun, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap=myblue)
# ax.gridlines()
ax.text(360-315,30,'g', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)


for drift in range(6): ax.plot([50+60*drift,50+60*drift],[90,40], transform=ccrs.PlateCarree(),c='lawngreen',linestyle='--')
for drift in range(4): ax.plot([35+90*drift,35+90*drift],[90,40], transform=ccrs.PlateCarree(),c='hotpink',linestyle='dotted')
cbar=fig2.colorbar(isfun,location='bottom',aspect=8,pad=0.01,shrink=0.85)
cbar.ax.set_xlabel('CH$_4$ fund [mW/m$^2$/sr]')

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text(((360-longit+180) % 360),40,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (40,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)



ax = plt.subplot(2,4,8,projection=ccrs.NorthPolarStereo(central_longitude=180))


ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
is323=ax.imshow(map_323, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),vmin=0,cmap='copper')
# ax.gridlines()
ax.text(360-315,30,'h', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)
cbar=fig2.colorbar(is323,location='bottom',aspect=8,pad=0.01,shrink=0.85)
cbar.ax.set_xlabel('F323N [mW/m$^2$/sr]')

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text(((360-longit+180) % 360),40,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (40,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)


# ax = plt.subplot(3,4,8,projection=ccrs.NorthPolarStereo(central_longitude=180))


# ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
# isdfun=ax.imshow(diff_methane_fun, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='seismic')
# ax.gridlines()
# for drift in range(6): ax.plot([50+60*drift,50+60*drift],[90,40], transform=ccrs.PlateCarree())
# cbar=fig2.colorbar(isdfun,location='bottom',aspect=8,pad=0.01,shrink=0.85)
# cbar.ax.set_xlabel('diff methane fund')
# ax.text(360-315,30,'d-ii', transform=ccrs.PlateCarree())

# ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
# ax.imshow(map_methane_hot, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap=myblue)
# ax.gridlines()


# ax = plt.subplot(3,4,9,projection=ccrs.NorthPolarStereo(central_longitude=180))
# ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
# isfun=ax.imshow(map_methane_hot, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap=myblue)
# ax.gridlines()

# for drift in range(6): ax.plot([50+60*drift,50+60*drift],[90,40], transform=ccrs.PlateCarree())
# cbar=fig2.colorbar(isfun,location='bottom',aspect=8,pad=0.01,shrink=0.85)
# cbar.ax.set_xlabel('methane hotband')
# ax.text(360-315,30,'e-i', transform=ccrs.PlateCarree())



# ax = plt.subplot(3,4,10,projection=ccrs.NorthPolarStereo(central_longitude=180))
# ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
# isdfun=ax.imshow(diff_methane_hot, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='seismic')
# ax.gridlines()
# for drift in range(6): ax.plot([50+60*drift,50+60*drift],[90,40], transform=ccrs.PlateCarree())
# cbar=fig2.colorbar(isdfun,location='bottom',aspect=8,pad=0.01,shrink=0.85)
# cbar.ax.set_xlabel('diff methane hotband')
# ax.text(360-315,30,'e-ii', transform=ccrs.PlateCarree())



# ax = plt.subplot(3,4,11,projection=ccrs.NorthPolarStereo(central_longitude=180))
# ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
# is323=ax.imshow(map_heat, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),vmin=0,cmap=myred)
# ax.gridlines()
# ax.text(180,50,'g')
# for drift in range(6): ax.plot([50+60*drift,50+60*drift],[90,40], transform=ccrs.PlateCarree())
# ax.text(360-315,30,'f-i', transform=ccrs.PlateCarree())


# cbar=fig2.colorbar(is323,location='bottom',aspect=8,pad=0.01,shrink=0.85)
# cbar.ax.set_xlabel('5.0-5.2$\mu$m')





# ax = plt.subplot(3,4,12,projection=ccrs.NorthPolarStereo(central_longitude=180))
# ax.set_extent([0, 360, 35, 90], ccrs.PlateCarree())
# isdfun=ax.imshow(diff_heat, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='seismic',vmin=-1,vmax=1)
# ax.gridlines()
# for drift in range(6): ax.plot([50+60*drift,50+60*drift],[90,40], transform=ccrs.PlateCarree())
# cbar=fig2.colorbar(isdfun,location='bottom',aspect=8,pad=0.01,shrink=0.85)
# cbar.ax.set_xlabel('diff 5.0-5.2$\mu$m')
# ax.text(360-315,30,'f-ii', transform=ccrs.PlateCarree())



fig2.savefig('asd_fig2b.pdf', dpi=300, bbox_inches='tight', facecolor='white', pad_inches=0) 

plt.show()



# %%


import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

fig2 = plt.figure(figsize=(12,8),dpi=300)


fig2.subplots_adjust(hspace=0,wspace=0.08)
bbox_props1= dict(boxstyle="circle,pad=0.15", fc="whitesmoke", ec="silver", lw=2)


ax = plt.subplot(1,1,1,projection=ccrs.NorthPolarStereo(central_longitude=180))
ax.set_extent([0, 360, 68, 90], ccrs.PlateCarree())
isfun=ax.imshow(map_methane_fun, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='blue_r',vmin=0.0)
# ax.gridlines()
ax.text(360-315,65,'a', transform=ccrs.PlateCarree(), ha="center", va="center", weight='bold',bbox=bbox_props1)
cbar=fig2.colorbar(isfun,location='bottom',aspect=8,pad=0.01,shrink=0.85)
# cbar.ax.set_xlabel('H$_3^+$ [mW/m$^2$/sr]')

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=1, color='k', alpha=0.5, linestyle='dotted')
gl.ylocator = mticker.FixedLocator([40,50,60,70,80])
for longit in range (0,355,60): ax.text((360-longit+180) % 360,70,str(longit)+'$^{\circ}$W', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)
for latit in range (70,82,10): ax.text(360-140,latit,str(latit)+'$^{\circ}$N', transform=ccrs.PlateCarree(), ha="center", va="center",size=8)


# ax = plt.subplot(242,projection=ccrs.NorthPolarStereo(central_longitude=180))
# # ax.set_extent([0, 360, 68, 90], ccrs.PlateCarree())
# # ax.imshow(diff_h3p, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='seismic',vmax=np.nanmax(diff_h3p),vmin=(-1)*np.nanmax(diff_h3p))
# ax.gridlines()
# # ax.text(310,50,'a')
# ax.plot([100,100],[0,50], transform=ccrs.PlateCarree())

